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Several mechanisms exist for generating a stochastic background of gravitational waves in the 
period following inflation. These mechanisms are generally classical in nature, with the gravitational 
waves being produced from inhomogeneities in the fields that populate the early universe and not 
quantum fluctuations. The resulting stochastic background could be accessible to next generation 
gravitational wave detectors. We develop a framework for computing such a background analytically 
and computationally. As an application of our framework, we consider the stochastic background 
of gravitational waves generated in a simple model of preheating. 



I. INTRODUCTION 

Prior to the time of recombination the universe is 
opaque to electromagnetic waves. The detection of a 
gravitational wave background generated before recom- 
bination would provide a unique observational window. 
Observations of this background would afford us rare and 
powerful probes of early universe physics and cosmology 
and could have profound implications. 

The detection and study of gravitational waves is cur- 
rently at the forefront of fundamental physics research. 
A world-wide network of gravitational wave detectors 
is poised to make the first detection. The Laser Inter- 
ferometer Gravitational- Wave Observatory (LIGO) and 
the Virgo interferometer have already achieved unprece- 
dented sensitivities at frequencies around 100 Hz and are 
currently undergoing upgrades. Future detectors such as 
Advanced LIGO and the Laser Interferometer Space An- 
tenna (LISA) promise to provide even more sensitivity 
and allow for observations at other frequencies. 

Inflation produces a stochastic background of gravi- 
tational waves through the amplification of primordial 
quantum fiuctuations [l|, 0]- Unfortunately this back- 
ground is too weak to be directly detected with existing 
instruments. It might, however, be observed via B-modes 
in the cosmic microwave background (CMB) 0,(3. In the 
period following inflation there are a number of mecha- 
nisms that could result in the production of an additional 
gravitational wave background; for example phase tran- 
To. [lli . as well as reheating and 

mIEM [13, m, [S US HH, m, m, 

24 . 125l l26j . These mechanisms are generally classical in 
nature: The gravitational waves are generated through 
inhomogeneities in the various matter flelds that popu- 
late the early universe. The gravitational waves produced 
by some of these mechanisms are close to being detectable 
by Advanced LIGO [13, [3 ; ^^'^ '^^^1 within the reach of 
LISA "27]. 



sitions d, 
preheatini 




'Electronic address: 
t Electronic address: 



larry ©gravity, phys . uwm. edu| 



slemens@gravlty.phys.uwm.edu| 



In many of these models the effects of the expan- 
sion of the universe cannot be neglected. The grav- 
itational waves can be produced on time and length 
scales comparable to the Hubble scale. In this paper we 
develop a framework for computing gravitational wave 
backgrounds in such situations and apply it to a simple 
well-studied model of preheating. 

The realization that preheating could lead to the effi- 
cient production of a gravitational wave background was 
first made by Khlebnikov and Tkachev [l^. They made 
an analytic estimate of the stochastic background pro- 
duced by preheating in a quartic inflation model. To esti- 
mate the stochastic background they used a flat space for- 
mula due to Weinberg [2^| for the energy in gravitational 
waves per unit solid angle. They found relatively large 
values of the background peaked at frequencies around 
/ ~ 10* Hz. A similar estimate by Garcia-Bellido p7j . 
showed that in hybrid inflation the frequency of the peak 
in the spectrum could be brought down to the 1 kHz 
range where Advanced LIGO might detect it. Almost 
a decade later, the problem was independently revisted 
by Garcia-Bellido and Figueroa [2l| and Easther and 
Lim [19]. The former [2l| studied hybrid inflation models 
using a flat space code that evolved the scalar fields as 
well as the metric perturbation. The latter [l^ evolved 
the scalar fields in a Friedmann-Robertson- Walker uni- 
verse using Felder and Tkachev's LatticeEasy [29] , con- 
sidered quartic as well as quadratic inflationary poten- 
tials, and used Weinberg's flat space formula. Easther, 
Giblin, and Lim [19, 23] later developed a new compu- 
tational strategy to include the effects of cosmological 
expansion. They coupled LatticeEasy to an integra- 
tor for the equation for the metric perturbation (using 
a mode decomposition) which they used to study hybrid 
inflation models. They found that the amplitude of the 
background is energy scale independent, and confirmed 
the analytic estimates of Garcia-Bellido [17] showing hy- 
brid inflation models might lead to a signal detectable by 
Advanced LIGO. Later, Garcia-Bellido, Figueroa, and 
Sastre ,20j re-examined the problem in an expanding 
spacetime using LatticeEasy and a configuration space 
integrator for the metric perturbation. In the mean- 
time Dufaux, Bergman, Felder, Kofman, and Uzan (2^ 
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studied the problem using Green's functions. They con- 
structed an approximate Green's function vahd for grav- 
itational waves generated well inside the horizon for any 
kind of cosmological expansion. The methods we develop 
in this paper are most similar to this work. 

In Section ini using a mode decomposition in the equa- 
tion for the evolution of the metric perturbation, we con- 
struct exact Green's functions for radiation- and matter- 
dominated expansions. Since our Green's functions are 
exact we faithfully recover the perturbations at all wave- 
lengths. We are, however, tied to particular types of ex- 
pansion (radiation- or matter-dominated). We use these 
Green's functions to construct an analytic expression for 
the energy emitted in gravitational waves as a function 
the Fourier transform of the stress energy tensor (in con- 
formal coordinates) analogous to Weinberg's flat space 
formula psj . This expression is useful when the evolution 
of the fields is known as well as for analytic estimates. 
We also use these Green's functions to write down the 
energy density in gravitational waves in a form suitable 
for lattice simulations. In Section IIIII we introduce the 
simple model of preheating that we use to validate our 
framework. This is followed by a prescription of how to 
translate the results of our simulations to values today. 
Finally we present and discuss our results. Our conclu- 
sions are foimd in Section ITVl 



II. THEORETICAL FRAMEWORK 

In this section we will lay out the basic framework 
for computing stochastic backgrounds given some mat- 
ter source (T)iiy). We will start by solving the perturbed 
Einstein equations for radiation- and matter-dominated 
expansion using exact Green's functions in conformal co- 
ordinates. In Appendix VK\ we provide these Green's func- 
tions in comoving coordinates. The metric perturbation 
determines the effective stress-energy tensor for gravita- 
tional waves, which, in turn, leads to an expression for 
the energy density in gravitational waves. We will de- 
velop this in both analytic and computational directions. 
On the analytic side we will provide a Weinberg-like for- 
mula for the energy in gravitational waves as a function of 
solid angle valid for expanding spacetimes. On the com- 
putational side, we will write down the same equations 
in a form useful for lattice simulations. 



A. Metric perturbations: Conformal Coordinates 

We will be primarily working in a spatially flat 
Friedmann-Robertson- Walker (FRW) metric in confor- 
mal coordinates: 



gauge (30| defined by 



a^{Tf){~dTf + dx^ + dy^ + dz^). 



(1) 



Our interest is in the first order metric perturbation. 
We choose to work in a spatial transverse-traceless 



TT 



0, 
0, 
0, 



(2) 
(3) 
(4) 



where is a spacetime index, i,j = 1,2,3 labels spatial 
indices and di = The physical metric is then 

ds^ = a^{ri)[-dTf + [5^j + hj^)dx'dx% (5) 

and the perturbed Einstein equations take on a particu- 
larly simple form: 



VTT , ^ a(??) ^rprp 



2i,TT 



(6) 



where the overdot denotes the derivative with respect to 
conformal time, is the three dimensional Laplacian 
in Euclidean space and the TT superscript denotes the 
transverse-traceless projection described below. 

For a generic perturbation of an arbitrary spacetime 
there is an inherent ambiguity in differentiating between 
the background and the perturbation. The ambiguity is 
essentially the gauge choice one makes when identifying 
points on the background spacetime with those on the 
physical (perturbed) spacetime [3l|. For a spatially flat 
FRW background, where the expansion is driven by a 
perfect fluid source with comoving velocity , the full 
( "background + perturbation" ) stress-energy tensor is 



(p + p)u^Uy + pg^^ + TTp 



(7) 



where p and p are the energy density and pressure, re- 
spectively, g^i, is the metric of the background (FRW) 
spacetime and tt^^ ~ tt^^ is the anisotropic stress-^. We 
consider the contribution of the anisotropic stress to be a 
perturbation of the otherwise homogeneous and isotropic 
background (with the expansion sourced by the first two 
terms in Eq. ([7])). The anisotropic stress tensor satisfies 
= 0, so it is both purely spatial and trace- 



less. In this situation it is clear that the only surviving 
term under the (spatial) transverse-traceless projection 
is ttJ^ , so this choice of gauge provides an unambiguous 
splitting of the (homogeneous) background and the per- 
turbation. For this paper we will be working exclusively 
in the TT gauge and will thus simply write our source 

term as T,^""". 

'■J . . 

Solving Eq. ([6]) is easiest if we first perform a spatial 
Fourier transform while leaving the dependence on con- 
formal time unchanged. Our convention for the Fourier 
transform is given by 



^) = 7^/ rf^rf'ke-^('"'-''--)/(;.,k). (8) 



Some authors write ttu 



in FRW backgrounds. 
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The perturbed Einstein equations then become 

(9) 

with = k • k. 

The transverse-traceless projection is easily described 
in momentum space. The TT part of some spatial tensor, 
Tij, is given by (S^l 

T^^^ik) = {P„„(k)F,„(k) - ipy(k)F„„(k)}T„„(k), 

(10) 



where 



kikj 



(11) 



77< = min(77, rjf), where 77 is the time at which the metric 
perturbation is evaluated and 77/ is the time the source 
turns off. 

It is interesting to note that both radiation- and 
matter-dominated expansion lead to relatively simple 
Green's functions. This is a result of the fact that we are 
working in momentum space. In configuration space the 
Green's function for radiation-dominated expansion only 
has support on the light-cone. The matter-dominated 
expansion Green's function contains a "tail" term that 
has support inside the past lightcone [1^. This severely 
complicates a configuration space calculation of the met- 
ric perturbation. 



B. Energy density in gravitational radiation 



In order to solve for the metric perturbation, we must 
know something about the evolution of the scale factor. 
We will restrict our attention to simple power law evolu- 
tion of the form 



0(77) = arj" 



(12) 



for some constants a and n. The solution in terms of a 
Green's function is given quite simply by 



167r 



drj' Lo 



In+l 



jn-l{uJT]')yn^l{uj'l]) 



The effective stress-energy tensor for gravitational ra- 
diation is given covariantly (in a transverse-traceless 
gauge) by 



/' 



(17) 



(13) 

where jniz) and yn{z) are spherical Bessel functions of 
the first and second kind, respectively. Our primary in- 
terest is in the particular cases 



0(77) = arj 



and 



a{ri) — arf 



(14) 



which correspond to radiation- and matter-dominated 
expansion, respectively. The general solution in Eq. (fT5|) 
then becomes 

/iJ.T(77,k) = / d77'77'sin[c.(77-^')]Tr(^''k), 

(15) 



where "f^v is the metric perturbation {g^^ = (^background _(_ 
7^1^) and the angle brackets denote a spatial average over 
several wavelengths and the covariant derivative is com- 
patible with the background metric. When working with 
this expression we must be careful to remember that ac- 
cording to Eq. (O 7JJ = a?hj!j^ . 

The quantity of interest for stochastic background 
computations is the energy density in gravitational 
waves, pgw. Though it is rarely stated explicitly, pgw 
is defined with respect to the proper time of a co-moving 
observer at rest, i.e., 

Pg„^t^i%s:(t,x), (18) 

where — (1,0,0,0) in the (background) metric 

ds^ = -di^ + a^{t)dx ■ dx. (19) 

Transforming to conformal coordinates: 



t^i'^TS-(t,x) 



'n 

for radiation-dominated expansion and 
167r /"'< 



/ir(^,k) 



' (w?7)3 



^77' 77'|(1 + Lu'^rirj') sin[cL;(77 — 77')] 
— Loiji — rj') cos[uj{ri — 77')] | 



_rfrf_^ 
32TTa^{r]) ^ 

i,3 



1 



327ra2(r/) 



5](/7TT(^,x)A5.T(,7,x)), (20) 



where rj^ — (1,0,0,0) in conformal coordinates. Per- 
forming the spatial averaging and using Parseval's theo- 
rem we obtain 



for matter-dominated expansion. The lower limit of in- 
tegration in both Eqs. and (fTB|) is determined by 
the time the source turns on, 77^. The upper limit is 



Pgw 



327ra(r;)2 V ^ 

id 



d-^k 



(21) 
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where V is the comoving volume over which the average is 
being performed. We will use this result in both analytic 
and computational contexts, which we develop separately 
in the next two subsections. 

It is important to note that not all modes of the pertur- 
bations we have been discussing are gravitational waves 
in the technical sense. Gravitational waves only exist as 
such when the typical wavelength of metric perturbations 
is much smaller than the characteristic length scale of 
the background spacetime, which in the present context 
is the Hubble radius. In other words, waves can only be 
observed in situations where a wave zone is well-defined. 
In the current situation, this only happens after the uni- 
verse has expanded significantly, e.g., today. Thus, in 
order to correctly apply the formalism developed in the 
remainder of this section, a (generally model dependent) 
procedure must be prescribed for "transferring" the value 
of at the time of creation to today's values (see Sec- 
tion nil Bl for an example) . 



C. Weinberg-like formula 

Perhaps the most often used tool for analytic calcu- 
lations of gravitational radiation is the the formula that 
appears in Weinberg [2^ for the energy in gravitational 
radiation, -Egw, as a function of solid angle^ 



dE, 



gw 



1,3 



dlO LO' 



\ ij 



(w,k)| 



(22) 



Though the formula is simple its usefulness is limited to 
applications in Minkowski space. It is straightforward 
to derive similar expressions for use in radiation- and 
matter-dominated (spatially flat) FRW spacetimes. 

Imagine that our source, T^"^ , is only active for some 
finite period of time before "turning off" . Then we are 
free to take the upper and lower limits of Eq. (fTS]) to 
infinity. At the same time, we can write the source as 
T^^{lu, k) using the inverse Fourier transform defined in 
Eq. ([S|). Interchanging the rj and lu integrals then leads 
to 



Lurj duj 



(u;,k), (23) 



where the details of the derivation can be found in Ap- 
pendix |BJ It follows directly that 

AS-^(,, k) ^ ^v/512^ ^^"("'^) -""J i-T^{-. k). 

(24) 



duj 



^ We have an extra factor of tt^ relative to the formula in Wein- 
berg's book [28j due to our Fourier transform convention [see our 
Eq. lH]. 



Using these expressions and setting a{rf) = arj in 
Eq. (ini), we have 

— [ duj [sin(a;77) — uirj cos(ti;?7)] ^ 
n ^ j J-oo 

2 



dEgv. 

dn 



o?rf 



(25) 



A similar expression holds for the case of matter- 
dominated expansion and any other simple power law 
evolution by using the Green's function in Eq. (|13p . 



D. Computational methods 

In many situations, such as the one we will consider 
in the next section, the source of gravitational radiation 
contains a nonlinear interaction term and we must resort 
to computer simulations for a result. The quantity of 
interest is given in conformal coordinates by 

, ,3 



dint 



327ra2(?7) V ^ 

1,3 



dn 



(r/,k) 



(26) 



This expression follows directly from Eq. (|2ip and is pro- 
portional to figw , the ratio of energy density in gravita- 
tional waves to the energy density required to close the 
universe. In practice, the integral over solid angle re- 
quires one to sum over the entire cubic lattice, which is 
an 0{N'^) operation, where N is the number of lattice 
points in each direction. 

In some situations, to reduce the computational com- 
plexity, a trick introduced in [2^ may be used. If we 
assume the stochastic background is isotropic, we can 
perform the integral analytically along any particular di- 
rection. That is, we can write 



lj2\h[;^{rj,cukp)\ (27) 



where kp is a unit vector in the direction we have cho- 
sen. In practice there are direction dependent statisti- 
cal fluctuations but these can be reduced by averaging 
Eq. (P7|) over several directions. There is a slight subtlety 
involved: For a cubic lattice with points, the maxi- 
mum length (in momentum space) of any component of 
k is proportional to N/2, so that in three dimensions the 
maximum possible length of k is Wmax oc \/3(N/2). How- 
ever, the length along any particular direction is generally 
shorter than this so the task is to try and maximize both 
'^max and the number of directions with that maximum 
length ijJmax- Following [2^1, we choose six face diagonals 
of the cube. Namely, the directions that run diagonally 
across the x-y, x-z, y-z, {—x)-y, (—x)-z and {—y)-z 
planes. This leads to a reasonable reduction in the effect 
of statistical fluctuations, with Wmax oc \/2{N/2). Fur- 
ther details of our computational methods are given in 
Appendix [Cl 
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III. EXAMPLE: PREHEATING 

In this section we will consider, as an application of our 
formalism, the stochastic background generated from a 
simple model of preheating after inflation. This is an 
active subject that has been studied with various tech- 
niques, which makes it ideal for validating our methods. 
We will start with a brief introduction to the physics of 
preheating. Then we will discuss how one uses the out- 
put of the simulation to determine the spectrum of the 
stochastic background today before providing a discus- 
sion of our results. 



0^ 








4 6 8 

Time (in units of the boxsize) 



A. Background 

At the end of inflation the universe is very cold and far 
from thermal equilibrium, which poses a challenge for the 
the hot big bang scenario. 

The earliest attempts at describing reheating after in- 
flation were centered on the idea that the oscillations of 
the inflaton about its minimum produce standard model 
particles that then interact and bring the universe into 
thermal equilibrium [H, IsJ, [H, [s^ . It was later discov- 
ered that this process is preceded by a stage of exponen- 
tial particle production which became known as preheat- 
ing [ilElla . 

We will consider a simple model of preheating in which 
the inflaton, (j), is coupled to some scalar field, Xi accord- 
ing to the Lagrangian 




V{cl^,x)], (28) 



where V{(f>, x) contains both the inflationary potential 
and the coupling between the two fields. In this work we 
will limit our consideration to 



mx) = ^0V + 



A 



(29) 



The coupling term between <j) and x in Eq. (|28p enters 
the Lagrangian as a spacetime dependent mass term for 
both fields. This non- linear term can produce an effect 
know as parametric resonance. Perhaps the most familiar 
example of parametric resonance is that of a child on a 
swing at the playground. The child, essentially a pendu- 
lum, learns fairly quickly if she swings her legs back and 
forth, effectively changing the length of the pendulum in 
a time-dependent manner, that she can swing higher and 
higher, having driven the system into resonance. 

A more precise description of the dynamics is given 
by which is partially summarized in the following. 
At the end of inflation, the x field starts off with zero 
amplitude and the oscillations of the inflaton about its 
minimum are described in terms of a conformally rescaled 
field If = a4> hy the equation 



(30) 



FIG. 1: Variances of (p and x fields, as computed by Lat- 
ticeEasy. Note that these units do not account for the ex- 
pansion of the universe. 



which has as its solution the Jacobi elliptic cosine func- 
tion 



1 



cn a: 



V2^ 



(31) 



where x — rjcjjoV^ and (/>o is the initial value of the in- 
flaton. The evolution of the /c-modes of the conformally 
rescaled Xk = axk is described by (the Lame equation) 



Xk 



A0g 



A 



cn X 



(32) 



This is essentially the harmonic oscillator with a time- 
dependent frequency given by the term in square brack- 
ets. The combination g^/X = q governs the strength of 
the time-dependent term in the frequency and is gener- 
ally called the resonance parameter. In our simulations 
we typically use q ^ 100, which results in complicated dy- 
namics. The coupling between x and the inflaton means 
that once the oscillations of the inflaton excite paramet- 
ric resonance in Xi it in turn brings the inflaton into a 
stage of parametric resonance, which continues until this 
process becomes inefficient. This is seen quite clearly in 
Fig. [1] which shows the variances of each of the fields as 
a function of time for one of our simulations (described 
in Section IIII Cp . Though the variances themselves are 
of no particular interest, they provide a simple way to 
visualize the time evolution of the full three dimensional 
dynamics. 



B. Values today 



As discussed in Section III Bl we need a prescription 
for translating metric perturbations in the early universe 
into gravitational radiation today. Our discussion of the 
transfer function, included for completeness, follows (2^ . 
I23j |. Gravitational energy density scales like a~'^ so that 



4 



(33) 
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where the and e subscripts denote today and the end 
of our simulations, respectively. Similarly, if entropy is 
conserved then radiation energy density scales like [3^ 



4 1/3 4 1/3 



(34) 



where go and gc are the number of relativistic species 
today and at the end time of our simulation (we take 
de/do = 100 for GUT scale inflation) and the energy 
density in radiation at the end of our simulations is, 
in the case of radiation-dominated expansion, just the 
total energy density at the end of the simulation (i.e., 
Prad.e — Ptot,e)- The quantity of interest is actually 



1 dp. 



Pclit d In uj 



(35) 



where Pcrit = ^Hq/Stt is the critical density required to 
close the universe and h is a dimensionless factor that 
accounts for the uncertainty in the value of the Hubble 
expansion rate today so that flg^jh^ is independent of 



this uncertainty. We can then use Eqs. 
write 



31 and (El to 




Aot,i 



dim 



(36) 



after inserting the appropriate values. Similarly, we want 
the physical frequency today, /o = wo/Stt — a;o/27rao, 
which is easily obtained using Eq. 



fo 



2nao 




1/4 
'^cPtot,c 



1/12 / V 1/4 

/ Prad,0 \ 

\ Piad,o / 
(4 X 10^° Hz). 



(37) 



The same expressions are found in [24I and [2S 



C. Results 

As stated previously, we use LatticeEasy for the evo- 
lution of the scalar fields. Conveniently, the two field 
model specified by Eq. pS)) is already setup in Lat- 
ticeEasy. For our comparison we used the same pa- 
rameters as many of the simulations published in |22l | 
and m. We took A = IQ-^-^ and g^ = 1.2 x lO-^^, 
which sets the resonance parameter q = 120. The initial 
values of the fields are 0o = 0.342Mp and xo = 0, where 
Mp is the Planck mass. This choice of parameters sets 
A(/)q/4 ~ (lO^^GeV)"*, which means inflation happens at 
the GUT scale. Additionally, LatticeEasy sets initial 




4 6 
Time (in units of the boxsize) 



FIG. 2: Evolution of the scale factor as a function of time, 
computed by LatticeEasy. Note that, to an extremely 
good approximation, the evolution is completely radiation- 
dominated. 



fluctuations of the fields as described in [29|. In confor- 
mal Planck units we took our box size to be L = 20, 
evolved the box for a time rjf = 240 using timesteps 
A77 — 0.01, and took TV = 256 points along each of the 
three axis (the typical size used in previous publications). 
We are able to run this simulation, including gravita- 
tional wave computations, in 18 hours on a four core 3 
GHz machine. The details of our implementation can be 
found in Appendix [Cl 

Figure [2] shows the self-consistent evolution of the 
scale factor as a function of time (the preferred unit for 
time in all of our plots is the box size of the simula- 
tion). The straight line indicates that the expansion is 
radiation-dominated throughout the course of the sim- 
ulation. More importantly, it validates our use of the 
radiation era Green's function. 

To compute the stochastic background spectrum we 
use Eq. (l36|) . The value of ptot.c is computed by Lat- 
ticeEasy. To estimate dpgw,o/rflnw we use Eq. (P7)) 
computed along six different directions as described at 
the end of Section fHDI Since the evolution is very well 
described by radiation-dominated expansion the values 
of hjj^ {ri,uj'kp) are updated at each time-step using 



167r 



drj' cos[ci;(77 — r]')] 

sin[(jj(?7 — 77')] I 

r] J 

xT;TT(^',c^kp) (38) 



which follows from Eq. (jT5|) . We do this integral using the 
rectangle method which appears to be sufficient. To com- 
pute Tj'^'^(?7',ti;kp) at each time-step we construct the six 
independent components of the spatial part of the stress 
energy tensor and Fourier transform them. Then we 
project out everything but the transverse-traceless part 
using Eqs. PTI)) and ITTI) for each of the six kp we have 
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Frequency [Hz] 

FIG. 3: The spectrum of the stochastic background today, 
Slgw/i^, computed along six directions on the lattice (thin col- 
ored lines), and the average (thick black line). 

previously chosen. The frequencies are given by Eq. ((37|) 
with We = 27r«-\/2/L, where i = 1 . . . N/2 (the factor of 
y/2 comes from our choice of kp along the face diagonals 
of the simulation box; see Section fll Dl for details). 

Fig. [3] shows the result of this procedure. We plot 
Q,fgv,h? as a function of frequency (both are the values 
today) . The thin colored lines are the results along each 
of the six directions in momentum space and the thick 
black line is the average. 

Next we present a comparison of our results with those 
obtained by other groups using different methods. The 
authors of [111, [1^, and [13] kindly provided us with 
the results of their simulations. The (solid) black curve 
is the result of our simulation (the same as in Fig. [3]). 
The (dash-dotted) blue curve is produced from the data 
published in ^23\^ . The (dashed) green curve shows the 
simulation results published in [23] for higher resolution 
simulation with N — 512 points in each direction. Fi- 
nally, the (dotted) red curve is from ^29] for a simulation 
with TV = 128^^ The results across all methods are in very 
good agreement. 

It is worth pointing out that the Green's function in- 
troduced in [22I ] is an approximation for general kinds 
of expansion, but for radiation-dominated expansion it 
is in fact exact. Looking at Eq. (12) of j2^] we see that 
it is the same as our Eq. (their metric perturbation 
is re-scaled by a factor of a{ii))- However, they work in 
the small wavelength approximation and perform an ad- 
ditional average over a single period of oscillation, which 
presumably accounts for the slight differences between 
the (dash-dotted) blue and (solid) black curves. 



^ These authors provided us with Qgw along each axis in addition 
to the six directions described in Section III Dl Only the latter 
was used in Fig. |4] 

* These authors provided us with Ogw/ptot,G and k at the end of 
the simulation. Fig. |4]displays the result of applying the transfer 
functions described in Section IlII BI to their data. 




Frequency [Hz] 

FIG. 4: Comparison of figw today as computed by the 
method described in this paper (solid black line) , as well as the 
curves published in [22] (dash-dotted blue line), [13] (dashed 
green line) and ^3] (dotted red line). Note the good agree- 
ment across methods. 



Although the results methods agree quite well with one 
another, there is a general issue common to all simula- 
tions: The fact that 77/ — 12L. This, by itself, is not a 
problem, but LatticeEasy uses periodic boundary con- 
ditions. This means that excitations of the scalar fields 
leave and re-enter the simulation box multiple times, set- 
ting up correlations in the field values not present in the 
early universe. The magnitude and importance of this 
effect is unclear. The problem can be avoided by setting 
the final time to (at most) half of the light-crossing time 
of the box, i.e., 77/ < L/2. Keeping the same final time, 
r]f — 240, we would set L — 480, which would require 
24^^ times as many points on the lattice to keep the cur- 
rent resolution. This, in turn would increase the memory 
required for the simulation by a similar factor. Alterna- 
tively, a glance at Fig. [T] suggests that gravitational wave 
production is essentially complete by 77 = 120 (this can 
be verified directly by monitoring the evolution in time 
of ilgw). A simulation with ?// = 120 and L = 240 would 
only require 12'^ times as many points to maintain the 
current resolution. We leave this investigation for future 
work. 

Finally, an obvious question that arises is what to 
do in situations where there is significant gravitational 
wave production during times when the expansion of the 
universe is not purely radiation- or matter-dominated. 
While each situation requires its own analysis, we can, 
generally speaking, combine our method with that of [2^ . 
That is, we can divide the calculation of flg„ into periods 
when the evolution is essentially a simple power law and 
use the exact Green's functions given by Eq. during 
those periods. In the periods in-between, where the evo- 
lution is either more complicated or changing rapidly, we 
can use the approximate expressions in [221 ]. 
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IV. CONCLUSIONS 

In this paper we have presented a method for comput- 
ing the stochastic background of gravitational waves due 
to (classical) sources in a spatially flat FRW background. 
Our method centers around "evolving" the metric pertur- 
bation using the exact Green's functions for radiation- 
and matter-dominated expansion. We developed these 
results into a Weinberg-like formula, applicable in ex- 
panding spacetimes, for use in situations where the dy- 
namics of the source is known or can be approximated. 
Additionally, we used the same results to present a com- 
putational framework applicable to scalar field sources 
with complicated dynamics. Because our framework re- 
lies on Green's functions for the evolution of the metric 
perturbation, it is robust and simple and we avoid the 
numerical difficulties associated with more complicated 
evolution schemes. Furthermore, we have shown that the 
results of using this framework to compute the stochas- 
tic background produced in the X(p'^ model of preheating 
agree quite well with those appearing in the literature. 

While the frequency ranges and sensitivities of the 
gravitational waves from the model of preheating stud- 
ied here make it unlikely to be detected in the near fu- 
ture, there remain several other options for the produc- 
tion of gravitational waves that could be observable soon. 
Aside from the sources discussed in Section |T] that may 
be observable to Advanced LIGO, phase transitions tak- 
ing place at the electroweak scale could lead to a signal 
detectable by LISA [s^. We plan to investigate this in 
more detail using the methods introduced here. 
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APPENDIX A: METRIC PERTURBATIONS IN 
COSMOLOGICAL TIME 

In this appendix we derive expressions for the metric 
perturbation in the commonly used coordinates 

ds^ = -dt^ + a^{t){6^J + hj^)dx'dx\ (Al) 

where hj^^ = hj^^{t,x) satisfies Eqs. ([IHl]). One could, 
in principle, perform a coordinate transformation on the 



expressions in conformal coordinates, but we find it more 
convenient to simply re-derive the results. In these coor- 
dinates the perturbed Einstein equations take the form 
(after performing a spatial Fourier transform): 



TT 



IGtt 



T 



TT 



(A2) 



Focusing once again on scale factor evolution of the form 
a{t) = ai", (A3) 
for constant alpha, leads to the Green's function solution 

1/2 



/■*< TT 



hl^it,k) = 16n 



i(l-3n)^/(l-n) 



2a2(n- 1) 



(A4) 



where Jpiz) and Yi3{z) are Bessel functions of the first 
and second kind, t< = min(f, tf) and 



1 - 3n 
2{n - 1) ' 

a{n — 1) ' 

ain — 1) 



(A5) 
(A6) 
(A7) 



Note that the factors in the denominators of the pre- 
ceeding expressions prevent them from being valid when 
n = I. The metric perturbation for a(t) — at is given in 
terms of 7 = -^/l — {uj/ay by 



dt' 



2ajt 



(A8) 



Particularly interesting cases of Eq. (jA4[) are radiation- 
and matter-dominated expansion which are characterized 
by 



2/3 



a(t) = at^/"^ and a(t) = at 



respectively. The corresponding metric perturbations are 
given by 



167r /■*< 
atot^^ A 



/iTT(i,k) = / di'sin 



2uj 



(^l/2_^,l/2) 



for radiation-dominated expansion and 



7^^(^',k), 
(A9) 
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16n 

9aLuH 



'U/3 



sm 



— (tl/3_t'l/3) 

a 



.^(il/3_i'l/3)cos 

a 



^(tl/3_i'l/3) 



(AlO) 



for matter-dominated expansion. 



APPENDIX B: DETAILS OF THE 
CALCULATION OF THE WEINBERG-LIKE 
FORMULA 

In this appendix we provide the details of the calcu- 
lation leading to Eq. ((23)) . from which the Weinberg- like 
formula in Eq. (I25p follows trivially. As mentioned pre- 
viously, we begin with Eq. (|15p and take the upper and 



lower limits of the integral to infinity while taking the in- 
verse Fourier transform of the source. At the same time 
we expand the sine into exponentials and make use of the 
identity 

r/ lim-i^e*'^"'. (Bl) 

Interchanging the order of the r\ and integrals then 
leads to Eq. (PH)) . More specifically 



167r 
^r] J 


' drj' T] 

— oo 


sin[a; 


Sni 




/>C30 


Lorj 


(-4) 


/ 

J —oo 


Stt 


Q roo 


duj' 



- — -^ r du;'i27Ty/^e^-^'S{u;'+u;-0 
dC 7-00 [ 



'<5(t^'-a;-C) T^"K,k) 



d 



Vl287r3 
ujr] 

Vl287r3 
Lorj 



_^^^sinM) a 



d 



^...,_j.TT(_^ + C,k) - e-'^-"-i;7(a; + C, k) 



e-|-T-(.,k)-e--A7;TT(^^k)) 



(B2) 



where we suppressed the limit symbol throughout, inte- 
grated by parts in the fifth line and took C ^ in the 
sixth line. From here it is straightforward to compute 
the derivative with respect to rj and use the expression 



for the energy density in gravitational radiation given 
by Eq. ([2T|) to obtain the Weinberg-like formula (quoted 
here for completeness): 



dEg^ 

In 



duj [sin(a;?7) — lot] cos{ujr])\ ' 



(B3) 



APPENDIX C: COMPUTATIONAL DETAILS lations on a four core 3 GHz machine. To take advantage 

Our code for computing the stochastic background in- 
terfaces directly with LatticeEasy. We ran our simu- 
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of th multi-processing capabilities available to us we used 
OpenMP |40i as implemented in gee version 4.2.3. The 
loops over the lattice necessary to evaluate the six inde- 
pendent components of Tij were split into four threads, 
each evaluated on one core. The most computationally 
expensive tasks we perform are the Fourier transforms 
of Tij [see Eqs. and (dH)]. There are six Fourier 
transforms to perform at each timestep. For this task we 
enlist the help of the Fastest Fourier Transform in the 



West (FFTW) Although FFTW has its own multi- 
processor capabilities, we find that we can get ^ 10% bet- 
ter performance using OpenMP for the task. By trial and 
error we found six threads, one for each Fourier trans- 
form, to be the fastest option. Note that we have not 
modified LatticeEasy. With these improvements the 
gravitational wave computation takes roughly the same 
amount of time as the field evolution. The simulation 
described in Section IIII CI takes 18 hours to complete. 
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